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ABSTRACT:  Molecular  mechanisms  of  the  action  of  antimicrobial  peptides  on  0  (IS  0.6  (IS  2.4  (IS 

bacterial  membranes  were  studied  by  large  scale  coarse-grained  simulations  of 
magainin  2— dipalmitoylphosphatidylcholine/palmitoyloleoylphosphatidylgly- 
cerol  (DPPC/POPG)  mixed  bilayer  systems  with  spatial  extents  up  to  0.1  /t m 
containing  up  to  1600  peptides.  Equilibrium  simulations  exhibit  disordered 
toroidal  pores  stabilized  by  peptides.  However,  when  a  layer  of  peptides  is  placed 
near  the  lipid  head  groups  on  one  side  of  the  bilayer  only,  their  incorporation  leads  to  a  spontaneous  buckling  of  the  bilayer.  This 
buckling  is  followed  by  the  formation  of  a  quasi-spherical  vesicular  bud  connected  to  the  bilayer  by  a  narrow  neck.  The  mean 
curvature  of  the  budding  region  is  consistent  with  what  is  expected  based  on  the  dependence  of  the  area  per  lipid  on  the  peptide-to- 
lipid  ratio  in  equilibrium  simulations.  Our  simulations  suggest  that  the  incorporation  of  antimicrobial  peptides  on  the  exterior 
surface  of  a  vesicle  or  a  bacterial  cell  leads  to  buckling  and  vesicle  budding,  presumably  accompanied  by  nucleations  of  giant  transient 
pores  of  sizes  that  are  much  larger  than  indicated  by  equilibrium  measurements  and  simulations. 


■  INTRODUCTION 

Antimicrobial  peptides  (AMPs)  are  found  in  many  organisms 
as  part  of  their  natural  defense  mechanisms  against  bacterial 
infection.  The  activity  of  AMPs  arises  from  their  interactions  with 
lipid  membranes  rather  than  proteins,  making  them  potential 
alternatives  to  traditional  antibiotics.  The  AMPs  are  usually  short 
and  highly  charged  cationic  a-helical  peptides  with  amphipathic 
character,  selectively  permeabilizing  bacterial  lipid  membranes  at 
sufficiently  high  concentrations  with  pore  formation.  The  exact 
nature  of  the  pore  formation  mechanisms,  however,  is  yet  to  be 
clarified  completely.  Among  the  models  proposed,1'2  the  carpet 
mechanism3,4  involves  peptides  adsorbed  on  the  membrane 
uniformly  up  to  a  critical  surface  concentration,  beyond  which 
a  global  instability  is  triggered.  The  barrel-stave  model5-7  has 
transmembrane  pores  formed  by  a  few  AMPs  lying  perpendicular 
to  the  membrane— water  interface  with  their  hydrophobic  faces 
stabilizing  the  pore  interior.  The  inner  surface  of  the  pore  is  lined 
up  instead  by  lipid  head  groups  in  the  toroidal  pore  model,8’9 
together  with  the  AMPs  lying  perpendicular  to  the  bulk  lipid 
surface.  In  the  sinking  raft  model,1  peptide  aggregates  on  the 
surface  sink  into  the  interior  maintaining  their  orientation 
perpendicular  to  the  interface,  eventually  reaching  the  other  side 
of  the  membrane.  Experimental  evidence  so  far  suggests  that  the 
action  of  different  AMPs  involves  varying  degrees  of  one  or 
multiple  mechanisms  among  these  possibilities. 

Magainin  2,  first  isolated  from  the  African  frog  Xenopus  laevis, 1 1 
is  one  of  the  most  well-characterized  AMPs.  In  its  interactions  with 
bacterial  membranes,  magainin  2  is  strongly  attracted  toward  the 
membrane— water  interface,  forming  an  a-helix  lying  parallel  to 
the  surface.1-’13  At  sufficiently  high  concentrations,  a  number  of 
these  helical  rods  are  believed  to  enter  the  membrane  interior,  in¬ 
ducing  pore  formation.14  Neutron  in-plane  scattering  experiments 


of  multilamellar  stacks8  and  vesicles15  suggested  that  these  inserted 
peptides  form  toroidal  rather  than  barrel-stave  pores,  and  the  mean 
diameter  of  these  pores  was  estimated  to  be  ~8  nm.15 

The  structure  of  equilibrated  pores,  however,  may  not  com¬ 
pletely  explain  the  mechanism  of  AMP  action  adequately,  since 
the  exact  transient  route  taken  by  systems  during  pore  formations 
likely  plays  important  roles  in  determining  the  fate  of  a  particular 
membrane  under  attack.  Indeed,  recent  studies  have  suggested 
that  kinetic  pathways  to  the  formation  of  such  pores  may  involve 
considerably  richer  structural  transformations  within  the  target 
membranes  than  suggested  by  equilibrated  pore  structures. 
Almeida  and  Pokorny  found  that  the  efflux  of  dyes  from  vesicles 
under  magainin  2  occurs  via  an  “all-or-none”  type  kinetic 
mechanism."  Giant  unilamellar  vesicles  (GUVs)  exhibit  expan¬ 
sions  of  surface  area  under  melittin  concentrations16  and  rapid 
leakages  of  fluorescent  dyes  under  magainin  2, 17,18  which  appears 
to  act  on  the  membrane  by  inducing  positive  spontaneous 
curvatures.18  When  the  composition  of  neutral-anionic  lipid 
mixtures  in  the  vesicles  is  varied,  the  minimum  peptide  concen¬ 
tration  required  for  pore  formation  increases  with  decreasing 
membrane  surface  charge  density.19  In  their  latest  study  using 
GUVs,  Tamba  et  al.20  found  that  the  time  course  of  leakage 
depends  on  the  size  of  fluorescent  probe  molecules.  With  small 
probes,  a  rapid  initial  leakage  was  followed  by  a  slower  regime.  In 
contrast,  larger  molecules  exhibited  a  much  smaller  initial  phase 
only  and  were  unable  to  leak  out  of  vesicles  completely.  The 
estimated  effective  cross  sectional  area  of  the  pores  correspond¬ 
ing  to  the  initial  leakage  was  ~20  nm. 
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These  results  together  imply  the  following  scenario  of 
magainin  2-induced  pore  formation:20  the  initial  strong  adsorption 
of  magainin  2  on  the  external  surface  of  the  membrane  causes  a 
large  imbalance  of  surface  tension  between  the  external  and 
internal  surfaces.  The  mechanical  instability  due  to  this  imbalance 
leads  to  the  creation  of  giant  pores  through  which  a  large  portion  of 
peptides  enter  the  vesicle  while  allowing  for  the  massive  effusion  of 
probes.  The  subsequent  balancing  of  the  concentration  of  peptides 
between  the  external  and  internal  surfaces  reduces  the  pore  size, 
which  eventually  reaches  the  level  observed  in  equilibrium  mea¬ 
surements  in  time  scales  up  to  minutes.8 

The  hypothesis  that  the  magainin  2-induced  pores  initially 
have  effective  diameters  much  larger  than  relevant  molecular 
length  scales  is  expected  to  have  significant  implications  to  its 
overall  action  in  vivo:  the  bacterial  cell  death  caused  by  the  pore 
formation  could  be  dominated  by  this  initial  transient  period,  in 
which  large  transient  forces  disrupt  the  membrane.  A  more  direct 
test  of  this  hypothesis  can  be  achieved  by  molecular  simulation 
techniques,  and  here  we  report  a  study  of  magainin  2-induced 
membrane  permeabilization  based  on  coarse-grained  (CG) 
molecular  dynamics  (MD)  simulations. 

Numerous  MD  simulations  have  previously  been  performed 
to  investigate  the  interaction  of  AMPs  with  membranes.  Appelt 
et  al.  have  studied  the  concentration  dependence  of  AMP— lipid 
interactions  using  all-atom  MD.21  Leontiadou  et  al.  performed 
MD  simulations  of  magainin-H2  in  dipalmitoylphosphatidylcho- 
line  (DPPC)  lipid  membranes,  reporting  the  formation  of 
disordered  toroidal  pores,  which  involved  up  to  four  peptides 
that  resulted  in  orientations  perpendicular  to  the  membrane 
surface.-2  This  study,  however,  employed  a  mutated  magainin 
that  interacts  strongly  with  neutral  bilayers.  Kandasamy  and 
Larson  have  investigated  the  effects  of  salt  concentrations  on 
magainin— bilayer  interactions  using  all-atom  MD.23  Dittmer 
et  al.  reported  a  combined  NMR— MD  study  of  alamethicin— 
lipid  interactions.24  Sengupta  et  al.  observed  significant  amounts 
of  disorder  in  toroidal  pore  structures  in  another  MD  study  of 
AMPs.25  Jean-Francois  et  al.26  also  investigated  the  roles  of 
electrostatic  effects  on  AMP  action. 

CG  simulations  greatly  expand  the  time  and  length  scales 
accessible  to  MD  studies  to  ranges  necessary  to  describe  large- 
scale  collective  effects  observed  in  experiments  with  high  peptide 
concentrations.  Iliya  and  Deserno  observed  via  solvent-free  CG 
simulations  that  concerted  interactions  between  multiple  AMPs 
can  stabilize  the  transmembrane  orientation.27  Other  AMPs  have 
also  been  investigated  in  a  number  of  recent  simulation  studies 
utilizing  CG  techniques,28  30  including  the  observation  of  the 
spontaneous  formation  of  barrel-stave  pores  in  LS3-DPPC 
systems  by  Gkeka  and  Sarkisov.31  Rzepiela  et  al.32  extended 
the  atomistic  simulation  study  of  magainin-H2/DPPC  system22 
with  combined  atomistic/CG  simulations  to  find  a  dynamic 
equilibrium  of  peptides  between  surface-aligned  and  transmem¬ 
brane  states  within  the  disordered  toroidal  pores. 

Major  advantages  of  CG  simulations  compared  to  atomistic 
studies  include  extensions  not  only  in  time  but  also  in  length 
scales.  However,  the  formation  and  evolution  of  mesoscale 
spatial  patterns  during  transient  pore  formations  induced  by 
AMP  action  have  not  been  examined  yet  using  CG  simulation  to 
our  knowledge.  In  this  paper,  we  report  a  CG  simulation  study  of 
AMP— bilayer  systems  in  length  scales  large  enough  to  examine 
the  transient  giant  pore  formation  suggested  by  recent  kinetic 
experiments.  We  used  the  empirical  CG  force  field  MARTINI 
developed  by  Marrink  and  co-workers,33  35  which  is  now  widely 


used  for  simulating  large-scale  collective  features  of  systems 
containing  membranes.  Gkeka  and  Sarkisov  have  recently  eval¬ 
uated  the  use  of  MARTINI  force  fields  in  characterizing  the 
interactions  between  amphiphilic  peptides  and  lipid  bilayers.36 
The  MARTINI  force  fields  for  lipids  in  particular  reproduce 
most  of  the  mesoscale  properties  of  monolayers  and  bilayers 
semiquantitatively,  including  the  gel— fluid  transition  of  DPPC 
molecules  (transition  temperature  295  ±  5  K  versus  the  experi¬ 
mental  315  K).37  Notably,  Baokina  et  al.38  have  elucidated  the 
molecular  mechanism  of  lipid  monolayer  collapse  under  lateral 
compression  using  the  MARTINI  force  fields,  revealing  the 
important  roles  played  by  the  buckling,  protrusion,  and  budding 
of  monolayer  folds  in  the  process. 

We  find  the  equilibrium  properties  of  magainin  2-containing 
bilayers  to  be  largely  consistent  with  previous  studies  that 
suggested  disordered  toroidal  pores  stabilized  by  peptides. 
However,  when  a  large  number  of  peptides  are  placed  only  on 
one  side  of  the  membrane,  we  observe  a  spontaneous  buckling  of 
the  bilayer  via  the  increase  of  positive  curvature  on  the  magainin 
2-bound  surface.  This  buckling  transition  is  akin  to  the  mono- 
layer  collapse  induced  by  lateral  compression  in  the  air— water 
interface.38  Our  results  provide  direct  structure-based  evidence 
for  the  important  roles  played  by  the  transient  giant  pore 
formation  during  AMP  action  as  suggested  by  recent  kinetic 
experiments.20 

■  METHODS 

Parameters  and  Simulation  Conditions.  Magainin  2 
(GIGKFLHSAKKFGKAFVGEIMNS)  and  mixed  DPPC/palmi- 
toyloleoylphosphatidylglycerol  (POPG)  bilayers  were  chosen  to 
build  a  model  of  an  AMP— anionic  lipid  membrane  system, 
simulated  using  the  MARTINI  force  field33,35  with  the  corre¬ 
sponding  parameters  for  DPPC  and  POPG.35,38  Na4  (or  Cl  ) 
ions  were  added  to  neutralize  each  system.  The  GROMACS 
software  package39  was  used  for  simulations.  The  molecules  in 
this  paper  refer  to  the  CG  MARTINI  versions  of  CG  groups.  The 
NMR  structure  of  magainin  2  in  micelles  by  Gesell  et  al.  (PDB  ID : 
2MAG)40  was  used  for  the  peptides.  The  DSSP  program41  was 
used  to  specify  the  secondary  structure  information  of  the  peptide, 
which  was  a-helical  for  18  residues.  Energy  minimizations  were 
performed  using  the  steepest  descent  method.  MD  simulations 
were  generally  performed  with  time  steps  of  0.03  ps  in  periodic 
boundary  conditions  with  electrostatic  options  recommended 
for  MARTINI’3,34  and  temperature  (310  K),  pressure 
(l  bar)  controls  using  Berendsen  thermostats.4"  The  pressure 
coupling  used  a  relaxation  time  of  3.0  ps  and  a  compressibility 
of  3  X  10-5  bar-1. 

System  Setup.  The  mixed  bilayer  containing  DPPC:POPG  in 
a  4:1  ratio  was  built  by  self-assembly  simulations.  Single  DPPC 
and  POPG  coordinates  were  replicated  to  96  and  24  molecules, 
respectively,  by  translations,  to  which  3200  water  and  24  Na4 
atoms  were  added  randomly.  The  resulting  system  was  energy 
minimized  and  simulated  for  27  ns,  which  produced  a  well- 
equilibrated  bilayer  configuration.  This  mixed  bilayer  was  repli¬ 
cated  together  with  a  peptide  placed  next  to  the  bilayer  interior. 
Alternate  inversions  of  z-direction  were  used  to  minimize  the  bias 
in  the  placement  of  peptides  toward  the  two  surfaces  of  the 
bilayer  to  build  a  system  with  one-fourth  of  the  compositions  of 
system  A  in  Table  1.  Equilibration  simulations  were  performed 
with  constant  pressure  setup  with  the  semi-isotropic  scheme  that 
allows  fluctuations  of  lateral  and  perpendicular  box  lengths 
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Table  1.  Configurations  of  Systems  Simulated 


number  of  molecules  dimension  (nm;  x  X  y  X  z) 


magainin  2 

DPPC 

POPG 

Na+/Cr 

CG  water11 

P/L 

simulation  length  ftts) 

initial 

final 

A” 

0 

6144 

1536 

1536/0 

204  800 

0 

0.6 

48.0  x  48.0  x  15.1 

77 

6144 

1536 

1305/0 

204  800 

0.0100 

0.6 

48.4  x  48.4  x  14.8 

162 

6144 

1536 

1050/0 

204  800 

0.0211 

0.6 

49.9  x  49.9  x  14.2 

230 

6144 

1536 

846/0 

204  800 

0.0299 

0.6 

50.3  x  50.3  x  13.0 

312 

6144 

1536 

600/0 

204  800 

0.0406 

0.6 

51.4  x  51.4  x  13.5 

384 

6144 

1536 

384/0 

204  800 

0.0500 

0.6 

51.5  x  51.5  x  13.5 

Bc 

384 

6144 

1536 

384/0 

554797 

0.0500 

2.4 

48  X  48  X  36 

37  x  37  x  60 

Cc 

384 

6144 

0 

0/1152 

554797 

0.0625 

2.1 

45  X  45  X  34 

34  x  34  x  66 

Dc 

1600 

24576 

6144 

1344/0 

2  301  023 

0.0521 

2.1 

95  X  95  x  36 

84  X  84  X  46 

Ed 

384 

5361 

1387 

235/0 

554797 

0.0569 

0.3 

37  X  37  x  60 

33  x  33  x  69 

“  A  CG  water  in  MARTINI  force  fields  corresponds  to  four  water  molecules.  h  Equilibrium  simulations  with  symmetric  distributions  of  peptides.  The  box 
dimension  refers  to  the  final  configurations  (Figure  l). 1  Transient  simulations  started  with  peptides  placed  on  one  side  of  the  bilayer  only  (Figure  4:  B; 
Figure  5:  C;  Figure  7:  D).  d  Configuration  of  the  system  containing  the  mini-vesicle  formed  by  removing  the  neck  in  the  final  configuration  of  B  (Figure  6). 


independently.  This  system  was  simulated  for  0.6  [is  and 
replicated  twice  in  the  x-  and  y-directions,  followed  by  another 
1.2  [is  MD  to  obtain  the  equilibrated  mixed  bilayer— magainin  2 
configuration.  The  resulting  system  A  with  the  peptide  to  lipid 
ratio  P/L  =  0.05  was  simulated  for  0.6  [is.  The  final  configuration 
was  used  to  generate  a  starting  configuration  of  a  system  with  a 
smaller  P/L  ratio  by  randomly  deleting  peptides  and  adjusting 
the  number  of  ions,  with  the  resulting  system  simulated  for 
0.6  [is.  This  procedure  was  repeated  up  to  P/L  =  0. 

Fluctuations  of  box  sizes  in  x-  and  y-directions  showed 
adequate  equilibrations  after  0.2  [is.  Equilibrium  properties  of 
magainin-embedded  bilayers  were  calculated  by  averaging  over 
0.2— 0.6  [is  for  each  simulation.  In  calculating  the  distribution  of 
magainin  2  for  Figure  2A,  the  peptide  center  was  defined  as  the 
center  of  mass  of  backbone  atoms.  The  height  was  defined  as  the 
difference  of  its  z-position  from  the  bilayer  center  at  the  location 
in  the  xyz-plane.  To  obtain  the  bilayer  center  height,  the  x,y-plane 
was  divided  into  grids  (20  in  each  direction,  or  ~2.6  nm  spacing), 
and  the  z-coordinates  of  phosphate  groups  within  each  grid  were 
averaged.  The  angle  of  peptide  orientation  in  Figure  2B  was 
defined  as  that  of  the  vector  connecting  the  two  backbone  atoms 
of  Glyl  and  Ser23  with  respect  to  the  z-axis  (modulo  90°).  The 
phosphate  group  distributions  in  Figure  3A  were  calculated  by 
similar  averages  of  phosphate  heights  with  respect  to  the  bilayer 
center  over  the  x,  y-plane  grids.  The  two  curves  for  P/L  =  0.05 
and  0  were  shifted  such  that  their  central  minima  would  be 
located  at  z  =  0. 

For  the  transient  membrane  buckling  simulations,  the  pure 
bilayer  ( P/L  =  0)  configuration  was  used  to  prepare  extended 
bilayers.  For  system  B,  a  slab  of  water  molecules  was  placed  on 
top  of  the  pure  bilayer  and  equilibrated  for  0.6  [is  prior  to  the 
placement  of  peptides.  To  build  system  D,  the  extended  bilayer 
used  for  B  was  replicated  twice  in  the  x-  and  y-directions  and 
equilibrated.  A  single  magainin  2  configuration  was  replicated 
into  x-  and  y-directions  with  spacings  chosen  such  that  the 
peptides  would  cover  the  membrane— water  interface.  Each 
peptide  was  then  translated  vertically  so  that  its  center  would 
be  positioned  1  nm  above  the  closest  phosphate  group  of  the 
bilayer,  producing  an  undulating  peptide  layer  that  uniformly 
covered  the  hydrophilic  region  near  the  solvent  interface.  Water 
molecules  closer  than  0.5  nm  to  peptide  atoms  were  then  deleted. 
The  systems  were  energy  minimized,  equilibrated  for  0.6  [is,  and 


simulated  for  periods  indicated  in  Table  1.  In  Figures  4,  5,  and  7,  a 
small  fraction  of  atoms  have  been  translated  based  on  the 
periodic  boundary  condition  such  that  the  bilayer  appears  as  a 
single  piece.  None  of  the  nonsolvent  atoms  were  closer  to  the 
periodic  images  of  others  along  the  z-direction  than  10  nm.  The 
pure  DPPC  bilayer— peptide  system  (system  C  in  Table  l)  was 
prepared  by  first  deleting  POPG  lipids  from  the  initial  config¬ 
uration  of  the  system  B  bilayer  without  peptides.  The  pure  DPPC 
bilayer  was  then  equilibrated  and  peptides  were  placed  similarly. 

For  the  simulation  of  vesicle  formation  from  the  disconnected 
bud,  the  lipid  molecules  in  the  neck  region  (between  10  and  20% 
of  the  box  length  in  the  z-direction  from  the  bottom)  were 
deleted,  and  the  number  of  ions  was  adjusted.  The  lipids  inside 
the  vesicle  were  classified  into  the  inside  and  outside  groups  by 
calculating  the  distance  of  phosphates  from  the  center  of  mass. 
The  radii  of  the  outside  and  inside  groups  were  obtained  by 
averaging  the  maximum  distance  from  each  phosphate  to  others 
in  the  same  group. 

■  RESULTS  AND  DISCUSSION 

Magainin  2  Forms  Disordered  Toroidal  Pores  in  Equilib¬ 
rium.  Equilibrium  model  bilayers  containing  DPPC:POPG  in 
4:1  ratio  and  magainin  2  peptides  were  built  and  simulated  in 
order  to  examine  the  stable  pore  morphologies  and  large-scale 
collective  effects  of  these  peptides  on  bilayer  properties.  The 
parameters,  lipid  composition,  and  temperature  (310  K)  are 
those  that  have  been  used  successfully  to  describe  lipid  mono- 
layer  collapse  by  Baoukina  et  al.38  In  building  these  systems, 
peptides  were  placed  in  the  interior  of  the  membrane  parallel  to 
the  lipid  tail  groups  before  equilibration,  such  that  they  would 
sample  relatively  rare  transmembrane  configurations  more  effi¬ 
ciently,  and  result  in  peptide  distributions  that  are  roughly  equal 
on  both  sides  of  the  membrane.  Figure  1  shows  an  example  of  the 
equilibrated  configurations  with  P/L  =  0.05,  which  corresponds 
to  a  surface  concentration  threshold  above  which  pore  forma¬ 
tions  on  GUVs  become  accelerated.19  Most  of  the  peptides  on 
the  surfaces  ended  up  in  clusters  containing  transmembrane 
pores,  which  appeared  to  perturb  the  bilayers  significantly 
(Figure  1A,B).  The  representative  pore  structure  shown  in 
Figure  1C  clearly  indicates  the  disordered  toroidal  morphology, 
where  peptides  lining  up  the  pore  wall  are  surrounded  by  a 
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Figure  1.  Equilibrated  bilayer  configuration  (DPPC:POPG  =  4:l) 
containing  magainin  2  at  the  peptide-to-lipid  ratio  of  P/L  =  0.05.  See 
Table  1  (system  A)  for  number  of  molecules  and  dimensions.  In  C,  a 
representative  pore  stabilized  by  the  peptides  is  shown  in  greater  detail. 
Color  schemes  are  as  follows:  P04  and  NC3  atoms  of  lipids  in  blue,  GL 
groups  of  POPG  in  red,  lipid  tail  carbon  atoms  in  green,  peptide  side 
chains  in  yellow,  and  peptide  backbone  groups  in  orange.  For  clarity,  CG 
water  groups  are  not  shown. 


Figure  3.  Modification  of  bilayer  properties  by  magainin  2.  (A)  Dis¬ 
tribution  of  phosphate  groups  within  the  membrane  for  peptide-free  and 
P/L  =  0.05  cases.  The  probability  distribution  is  normalized  such  that  the 
total  area  under  the  curve  equals  1.  (B)  Dependence  of  area  per  lipid  on 
P/L.  The  dashed  line  is  a  linear  fit  (area/nm2  =  0.598  +  2.081  P/L). 


A  B 


Figure  2.  Equilibrium  properties  of  the  pores  formed  in  magainin  2-im- 
bedded  bilayers.  A  and  B  show  the  distributions  of  peptide  center  and 
orientation  within  the  bilayer,  respectively,  at  P/L  =  0.05  (see  Figure  l). 


substantial  amount  of  lipid  head  groups  drawn  into  the  bilayer 
interior.  The  distribution  of  peptide  centers  of  mass  (Figure  2A) 
suggests  that  at  peptide  concentrations  relevant  to  membrane 
disruption,  magainin  2  are  distributed  widely  over  the  pore 
interior  including  its  center.  On  the  other  hand,  their  orientation 
with  respect  to  the  bilayer  surface  (Figure  2B)  exhibits  a  bimodal 
distribution  with  two  possible  orientations,  one  parallel  to  the 
surface  (~90°),  and  another  near  ~30°.  The  distribution  near 
the  orientation  parallel  to  the  surface  is  sharp  but  smaller  in 
fraction,  consistent  with  the  expectation  that  they  would  appear 
within  the  tail  parts  of  the  positional  distribution  of  Figure  2A. 

The  incorporation  of  peptides  into  the  bilayer  in  the  form  of 
these  toroidal  pores  leads  to  modifications  of  its  overall  equilibrium 
properties.  The  distribution  of  lipid  head  parts  (phosphate  groups) 
across  the  bilayer  (Figure  3A)  indicates  that  as  magainin  2 
molecules  bind  and  form  pores,  the  bilayer  becomes  thinner  while 
increasing  the  density  of  head  groups  at  the  bilayer  center, 
consistent  with  the  toroidal  nature  of  pore  morphology  suggested 
in  Figure  1C.  The  thinning  ofbilayer  thickness  upon  AMP  addition 
has  widely  been  observed  experimentally,16’43  including  magainin 
2-containing  bilayers.44  The  main  cause  of  such  thinning  effects  is 
the  expansion  of  the  total  surface  area  of  bilayers  containing  AMPs, 
which  we  also  directly  observe  in  Figure  3B,  where  the  area  per 
lipid  is  seen  to  increase  linearly  with  increasing  P/L. 

Our  large-scale  equilibrium  simulation  results  for  P/L  =  0.05 
are  consistent  with  previous  experimental  and  simulation  studies 
that  provided  evidence  for  the  formation  of  disordered  toroidal 
pores  at  high  concentrations  of  magainin  2  in  bilayers  with  net 
negative  surface  charge  densities.  Ludtke  et  al.  performed 
neutron  in-plane  scattering  experiments8  to  show  that  the  pores 
formed  by  magainins  are  much  larger  in  size  than  can  be 


explained  by  the  barrel-stave,  proposing  the  toroidal  pore  model. 
In  atomistic  MD  simulations  by  Leontiadou  et  al22  of  magainin 
H2— DPPC  systems  containing  up  to  four  peptides,  spontaneous 
formations  of  such  disordered  toroidal  pores  were  observed 
within  0.25  /is.  The  equilibrium  pore  structure  we  observe 
(Figure  IC)  is  largely  consistent  with  the  previously  noted 
atomistic  as  well  as  the  combined  atomistic-CG  results.32  The 
latter  study,  in  particular,  showed  a  similar  increase  in  phosphate 
group  density  at  the  bilayer  center  as  well  as  the  thinning  of 
bilayer  thickness  on  magainin-H2  addition  to  the  neutral  DPPC 
bilayer.  Our  results  greatly  extends  the  length  scale  of  the  system 
under  consideration  (from  ~4  peptides  to  ~400)  and  show  that 
the  equilibrium  properties  observed  in  CG  simulations  using 
MARTINI  force  fields  remain  applicable  with  more  experimen¬ 
tally  relevant  magainin  2— anionic  lipid  bilayer  systems.  Never¬ 
theless,  the  kinetics  of  pore  formation  we  deduce  from  our  study 
(see  below)  is  quite  distinct  from  what  was  suggested  in  the 
simulation  of  the  modified  magainin-H2  and  neutral  bilayer. 

Asymmetric  Magainin  2  Incorporation  Leads  to  Buckling. 
When  a  bacterial  cell  is  exposed  to  high  concentrations  of  AMPs, 
the  peptides  will  initially  have  to  bind  to  the  outer  layers  of  the 
membrane  only.  This  asymmetry  is  expected  to  cause  significant 
differences  in  the  kinetic  behavior  ofbilayer  morphology  changes 
compared  to  what  is  observed  in  equilibrium  conditions.  In  order 
to  simulate  these  kinetic  changes  and  infer  their  implications  to  the 
mechanism  of  AMP  action,  we  prepared  a  system  of  mixed  bilayer 
with  a  large  vertical  extension  and  placed  a  layer  of  magainin 
2  molecules  on  top  of  one  side  of  the  lipid— water  interface 
(Figure  4).  The  peptides  quickly  became  incorporated  into  the 
headgroup  regions  with  a  moderate  degree  of  aggregation.  How¬ 
ever,  the  bilayer  soon  developed  a  positive  curvature,  with  the 
curved  regions  of  the  bilayer  eventually  coalescing  into  a  vesicular 
bud  in  ~0.6  /is.  The  budding  protrusion  was  coated  with  peptides 
on  the  outer  surface  and  remained  connected  to  the  bilayer  with  a 
neck.  The  growth  of  the  protrusion  was  essentially  complete  within 
~1  /is,  and  the  system  remained  stable  for  another  ~1  /is.  During 
this  buckling/budding  process,  the  dimension  of  the  system 
changed  with  a  large  decrease  and  increase  in  sizes  in  the  x-,  y-, 
and  z-directions,  respectively  (Table  l). 

The  peptides  remained  strongly  associated  with  membrane 
head  groups  in  all  of  our  simulations.  This  strong  attractive 
interaction  allows  our  simulation  setup  to  maintain  the  asym¬ 
metry  of  peptide  concentrations  across  the  bilayer  even  under 
severe  distortions.  Alternative  setups  such  as  placing  two  bilayers 
within  the  unit  cell  under  periodic  boundary  conditions  often 
used  for  simulating  ion  transport  are  therefore  not  necessary  for 
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Figure  4.  Vesicular  budding  from  the  bilayer  (DPPC:POPG=4:l)  with 
magainin  2  (yellow)  adsorbed  on  one  side  (System  B  in  Table  l).  For 
each  time  snapshot,  the  bottom  figure  shows  the  cross-sectional  view. 


our  system.  However,  we  did  find  that  large  scale  changes  in  the 
aspect  ratio  accompanying  the  buckling  process  may  produce 
configurations  in  which  periodic  images  of  one  monolayer 
approach  the  other,  which  can  allow  peptide  redistributions 
across  the  bilayer.  We  therefore  chose  the  initial  configurations 
for  the  transient  simulations  to  be  sufficiently  large  in  order  to 
avoid  such  artifacts. 

The  spontaneous  buckling  observed  with  asymmetric  bind¬ 
ing  of  peptides  has  previously  been  postulated  based  on  many 
kinetic  experiments  of  AMP  action  on  negatively  charged 
bilayers.16’17'19,20  Tamba  and  Yamazaki  observed  that  GUVs 
with  nonspherical  (such  as  prolate)  shapes  develop  narrowing 
necks  and  evolve  into  dumbbell  and  necklace-shaped  objects 
when  exposed  to  increasing  concentrations  of  magainin  2  on 
their  exterior  surfaces.17  Such  shape  changes  can  be  explained 
by  the  increase  in  effective  surface  area  of  the  outer  membrane 
due  to  the  binding  of  peptides.  In  equilibrium  experiments44 
and  simulations  (Figure  3),  the  expansion  occurs  symmetrically 
on  both  sides  of  the  bilayer,  which  preserves  its  planar  average 
geometry.  In  contrast,  when  the  area  increase  is  induced  only  on 
one  interface  of  the  bilayer  as  in  a  nonspherical  GUV  placed 
inside  a  magainin  2  solution,  the  difference  in  area  across  the 
bilayer  must  be  accommodated  by  shape  changes  under  the 
constraint  of  constant  enclosed  volume.  These  morphological 
changes  can  be  described  effectively  by  the  continuum  elastic 
descriptions  of  vesicles  such  as  the  area-difference  elasticity 
model,45  which  predicts  budding  transitions  via  the  growth  of 
protrusions,  narrowing  of  necks,  and  the  eventual  fission  of 
the  daughter  cells  as  has  indeed  been  observed  in  GUV 
experiments.17  Lee  et  al.  also  have  directly  noted  morphology 
changes  of  GUVs  in  melittin  solution,  where  protrusions  developed 
into  micropipets  holding  the  vesicles.16 

Our  transient  buckling  simulations  confirm  the  driving  force  of 
such  morphology  changes  (unbalanced  expansion  due  to  asym¬ 
metric  peptide  adsorption)  by  connecting  molecular  level  fea¬ 
tures  to  the  mesoscale  phenomenon.  The  periodic  boundary 
condition  imposed  in  the  simulation  replaces  the  constant 
volume  constraint  of  GUVs  in  experiments.  The  growth  of  the 
protrusion  into  a  vesicular  bud  connected  to  the  base  bilayer  that 
roughly  remains  planar  can  be  interpreted  as  a  result  of  the 
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Figure  S.  Vesicular  budding  from  the  pure  DPPC  bilayer  with  magainin 
2  placed  on  one  side  (system  C  in  Table  l). 


optimization  of  the  free  energy  under  this  boundary  condition 
constraint. 

To  examine  the  extent  to  which  this  buckling  process  depends 
on  lipid  compositions,  we  have  performed  a  similar  simulation 
with  a  bilayer  consisting  of  DPPC  only  (Figure  5).  The  mor¬ 
phological  shape  change  and  its  time  scale  are  similar  to  those 
observed  for  the  charged  bilayer  in  Figure  4.  This  result  suggests 
that  the  global  net  charge  of  the  bilayer  does  not  play  important 
roles  for  the  buckling  process  simulated,  where  a  planar  layer  of 
peptides  are  placed  near  the  headgroup  of  a  monolayer,  mimick¬ 
ing  high  initial  surface  concentrations  of  peptides.  In  GUV 
experiments,  decreasing  the  surface  charge  density  of  the 
bilayer  (or  increasing  the  salt  concentration)  leads  to  increases 
both  in  the  magainin  2  concentration  in  the  buffer  required  for 
the  pore  formation  as  well  as  its  time  scale.19  The  fact  that  the 
artificial  placement  of  peptides  near  the  monolayer  headgroup 
leads  to  buckling  irrespective  of  lipid  composition  implies  that 
the  selectivity  of  magainin-induced  pore-formation  for  charged 
membranes  is  largely  based  on  the  net  long-range  attraction 
between  individual  peptides  and  the  bilayer  rather  than  their 
collective  action  near  contact.  For  neutral  or  weakly  charged 
bilayers,  the  gradual  adsorption  and  accumulation  of  peptides  on 
the  outer  monolayer  as  modeled  by  the  initial  configurations  of 
Figures  4  and  5  would  either  be  impossible  or  require  higher 
peptide  concentrations  and  time  scale. 

Experiments  have  also  indicated  that  the  key  parameter  that 
determines  the  rate  of  pore  formation  is  the  mole  fraction  Xj,  of 
peptides  within  the  bilayer  rather  than  the  ratio  of  bulk  magainin 
concentration  to  lipid  concentration,  with  the  critical  ~ 
0.05. 19  However,  since  in  our  simulations  practically  all  peptides 
remain  tightly  incorporated  into  the  bilayer,  the  P/L  value  we 
used  should  be  interpreted  as  a  measure  of  surface  rather  than 
bulk  concentration  (Xj,  =  0.0476  for  Figures  4  and  7).  Our 
simulated  system  therefore  roughly  corresponds  to  the  threshold 
regime  in  GUV  experiments  where  the  effective  peptide  con¬ 
centration  in  the  bilayer  begins  to  cause  a  large  increase  in  the 
leakage  rate. 

A  Spherical  Mini-Vesicle  Forms  When  the  Neck  Is  Re¬ 
moved.  The  stabilization  of  the  neck  thickness  to  a  relatively 
large  value  reflects  the  fact  that  the  fission  of  vesicles  across  the 
point  of  singularity  in  continuum  descriptions45  involves  a  high 
free  energy  barrier  associated  with  a  simultaneous  reversal  of 
orientations  of  a  large  number  of  lipid  molecules  inside  the 
bilayer.  This  activation  process  usually  is  catalyzed  by  proteins  in 
cellular  in  vivo  conditions.46  These  proteins  (clathrins)  assemble 
themselves  into  coats  on  the  surface  of  the  bilayer,  inducing 
budding  and  fission.  Kozlovsky  and  Kozlov  studied  the  ener¬ 
getics  of  this  fission  process  using  elastic  models  and  estimated 
the  barrier  to  be  on  the  order  of  ~  100  kBT  47  The  CG  simulation 
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Figure  6.  Spontaneous  stabilization  of  a  mini-vesicle  from  a  discon¬ 
nected  bud.  The  initial  configuration  (A,  B  each  showing  the  side  and 
bottom  views)  was  generated  by  deleting  lipid  molecules  in  the  neck 
region  from  the  configuration  at  2.4  [is  of  Figure  4  (shown  here  rotated 
clockwise  by  90°).  C  and  D  show  the  configuration  after  0.3  [is.  See 
Table  1  for  the  system  compositions  and  dimensions  (System  E). 


result  shown  in  Figures  4  and  5  suggests  that  AMPs  including 
magainin  2  may  act  in  ways  similar  to  these  fission  proteins  by 
bending  the  bilayer  into  small  buds.  The  actual  fission  of  these 
buds,  however,  would  have  to  be  activated  by  rare  but  non- 
negligible  thermal  fluctuations. 

To  examine  whether  such  mini- vesicles  with  diameters  of  tens 
of  nanometers  can  indeed  form  as  a  result  of  AMP  action  when 
the  fission  free  energy  barrier  is  overcome,  we  removed  lipid 
molecules  in  the  neck  region  of  the  vesicular  bud  (Figure  6).  Both 
the  disconnected  bud  and  the  base  membrane  rapidly  healed 
their  openings,  producing  a  spherical  mini-vesicle  coated  with 
peptides  and  a  bilayer  containing  the  remaining  fraction  of 
magainin  2.  The  resulting  vesicle  contained  239  peptides,  and 
1280  and  1808  lipids  in  the  internal  and  external  surfaces, 
respectively.  The  radii  of  the  internal  and  external  monolayers 
were  rin  =  7.56  nm  and  rout  =  1 1.40  nm,  respectively.  The  area  per 
lipid  is  therefore  4 71  ( r m2  +  rout2)/(l280  +  1808)  =  0.761  nm2, 
and  P/L  =  0.077.  The  equilibrium  relationship  shown  in 
Figure  3B,  on  the  other  hand,  predicts  0.758  nm"  for  the  same 
P/L,  in  good  agreement  with  the  actual  geometry.  For  compar¬ 
ison,  the  initial  configuration  of  the  flat  bilayer  in  Figure  4  had  an 
area  per  lipid  of  ~0.59  nm2.  The  fact  that  the  spherical  vesicle 
derived  from  our  CG  simulations  started  from  a  flat  bilayer  has 
geometries  expected  from  equilibrium  properties  further  supports 
the  conclusion  that  the  spontaneous  buckling  and  subsequent 
morphological  evolution  are  driven  by  the  energetic  penalty 
associated  with  the  surface  tension  imbalance  across  the  bilayer. 

The  Size  of  Vesicular  Bud  Is  Independent  of  System  Sizes. 
In  order  to  test  how  the  simulated  system  size  affects  this 
protrusion  process,  we  examined  the  time  evolution  of  a  system 
that  is  twice  the  size  in  the  x-  and  y-directions  at  the  same  P/L 
value  (Figure  7).  A  buckling  process  similar  to  that  in  Figure  4 
was  observed  over  the  course  of  ~1  [is,  which  again  produced  a 
vesicular  bud  after  ~2  [is  with  a  radius  of  ~10  nm.  The  fact  that 
doubling  the  size  of  the  overall  simulation  cell  with  the  same  P/L 
value  leaves  the  size  of  the  vesicular  bud  unchanged  suggests  that 
the  periodic  boundary  condition  imposed  on  simulations  does 
not  substantially  alter  the  physical  mechanism  behind  the  bud¬ 
ding  process. 

This  larger  system  had  the  spatial  extension  approaching 
0.1  [im  for  the  size  of  bilayer  and  contained  1600  magainin  2 
peptides,  making  it  roughly  2  orders  of  magnitude  larger  than  any 
previous  simulation  studies  of  AMP— bilayer  systems.  This 
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Figure  7.  Vesicular  budding  from  a  peptide-containing  bilayer  system 
twice  as  large  as  that  in  Figure  4.  See  Table  1  for  system  composition  and 
dimensions  (system  D). 


spatial  extent  nevertheless  still  covers  less  than  ~10%  of  a  single 
GUV  or  a  bacterial  cell  that  is  on  the  order  of  ~1  [im  in  size. 
Flowever,  for  spherical  GUVs  for  which  global  morphological 
changes  into  shapes  such  as  dumbbells  or  necklaces  are  not 
directly  accessible,  it  is  expected  that  the  poration  process  will 
occur  locally  on  the  surface  while  the  overall  spherical  morphol¬ 
ogy  of  the  vesicle  is  still  maintained:  experimental  studies  of 
GUVs  of  ~10  [im  diameter20  revealed  no  discernible  changes  to 
the  vesicle  shape  even  after  a  profuse  leaking  has  already  been 
initiated.  Our  simulation  results  suggest  a  possible  mechanism  in 
which  such  transient  pores  may  initially  form  on  the  surface:  the 
growing  imbalance  in  surface  tensions  between  the  exterior  and 
interior  surfaces  from  peptide  binding  begins  to  instigate  buck¬ 
ling  of  bilayers.  The  volume  of  the  vesicle  and  the  average  surface 
area  of  the  enclosing  sphere,  however,  are  constrained,  which  do 
not  allow  the  contraction  of  the  sort  evident  in  simulations 
shown  in  Figures  4,  5,  and  7  where  box  sizes  in  the  x-  and 
y-directions  shrink  as  the  vesicular  protrusion  develops.  The 
buckling  and  protrusion,  nevertheless,  can  occur  if  pores  of 
similar  sizes  are  formed  in  neighboring  regions  within  the  bilayer. 
Such  pore  formations  in  bilayers  require  crossing  the  nucleation 
barrier  beyond  which  the  line  tension  is  overcome  by  the 
membrane  tension.48  The  AMPs  on  the  surface  can  lower  the 
tension  and  reduce  the  critical  radius.  A  rare  thermal  activation 
would  still  be  necessary  to  cross  the  barrier,  which  is  consistent 
with  the  experimental  observation  that  the  initiation  of  leakage 
on  single  GUVs  under  high  concentrations  of  magainin  2  is 
highly  stochastic,  with  waiting  times  on  the  order  of  ~  100  s.20 

We  have  also  simulated  the  same  system  under  constant  area 
(box  sizes  in  the  x-  and  y-directions  fixed)  conditions  to  see  to  what 
extent  such  a  nucleation  of  pores  might  be  accessible  via  uncon¬ 
strained  simulations.  A  moderately  perturbed  bilayer  shape  (similar 
to  the  0.3  [is  configuration  in  Figure  7 )  persisted  up  to  ~  1 .2  [is  with 
litde  changes  in  its  degree  of  buckling,  consistent  with  a  high  free 
energy  barrier  of  pore  nucleation.  We  note,  however,  that  many 
other  factors  beyond  the  reach  of  CG  simulations  may  contribute 
significantly  to  the  pore  formation  mechanism:  for  instance,  during 
the  long  waiting  time  periods  (~100  s)  prior  to  the  activated  pore 
formation  event,  appreciable  amounts  of  water  may  diffuse  across 
the  membrane,  allowing  for  changes  in  GUV  volume  and  surface 
area.  Simulation  studies  using  different  schemes  of  boundary 
conditions  designed  to  probe  such  effects  may  offer  further  insight 
into  the  precise  mechanism  of  pore  formation. 
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Nevertheless,  it  is  particularly  notable  that  the  experimen¬ 
tally  estimated  mean  radius  of  the  giant  transient  pores  that 
form  initially  (~20  nm)20  is  similar  to  the  length  scale  of 
vesicular  protrusions  we  observe  in  Figures  4  and  7.  In  the 
GUV  experiments,  the  bulk  of  leakage  occurs  after  the  initial 
giant  pore  formation,  followed  by  the  gradual  decrease  of  the 
pore  radius  to  levels  observed  in  equilibrium  measurements. 
The  two  very  different  classes  of  bilayer  disruption  we  show  in 
Figure  1  and  in  Figures  4,  5,  and  7,  one  in  equilibrium  where  the 
peptides  are  symmetrically  distributed  within  the  toroidal 
pores,  and  the  other  where  they  cause  violent  buckling/ 
budding  (accompanied  by  giant  pore  formations),  provide 
strong  supporting  evidence  for  the  hypotheses  put  forward 
in  recent  kinetic  experiments. 16,19,20 

■  CONCLUSIONS 

Although  equilibrium  structures  of  pores  formed  by  AMPs 
in  bilayers  have  been  studied  widely  both  in  experiments  and 
simulations,  the  nature  of  kinetic  routes  a  vesicle  or  a  bacterial 
cell  experience  when  under  attack  by  such  peptides  has 
received  relatively  less  attention.  Recent  experiments  using 
single  GUVs  have  suggested  that  the  initial  pore  formation 
responsible  for  the  bulk  of  leakage  involves  length  scales  much 
larger  than  those  indicated  by  the  equilibrium  structures.  In 
our  study,  we  greatly  extended  the  length  scale  of  simulations 
from  levels  of  previous  works  on  AMPs  such  that  the  relevant 
mesoscale  collective  effects  can  be  observed  directly.  The 
results  obtained  provide  direct  structural  and  kinetic  evidence 
supporting  the  roles  played  by  such  nonequilibrium  effects  in 
AMP  action:  the  imbalance  of  surface  concentrations  across 
the  bilayer  provides  the  driving  force  for  large  scale  structural 
instabilities,  which  may  lead  to  buckling,  budding,  and  nuclea- 
tion  of  giant  pores. 
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